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Two-dimensional Josephson junction arrays at zero temperature are investigated numerically 
within the resistively shunted junction (RSJ) model and the time-dependent Ginzburg-Landau 
(TDGL) model with global conservation of current implemented through the fluctuating twist 
boundary condition (FTBC). Fractional giant Shapiro steps are found for both the RSJ and TDGL 
cases. This implies that the local current conservation, on which the RSJ model is based, can be re- 
laxed to the TDGL dynamics with only global current conservation, without changing the sequence 
of Shapiro steps. However, when the maximum widths of the steps are compared for the two models 
some qualitative differences are found at higher frequencies. The critical current is also calculated 
and comparisons with earlier results are made. It is found that the FTBC is a more adequate 
boundary condition than the conventional uniform current injection method because it minimizes 
the influence of the boundary. 

PACS numbers: 74.50.+r, 74.40.+k, 74.60.Ge, 74.60.Jg 



Two-dimensional (2D) Josephson junction arrays 
(JJA's) have been the subject of much current inter- 
est because of their phase transitions and noncquilibrium 
transport properties!!] In the presence of an external uni- 
form magnetic field, the frustration /, defined by the 
number of flux quanta per plaquette, plays an important 
role reflected, e.g., in the value of the critical current.B 
Furthermore, when an L x L square array with / = p/q 
(p, q are integers) is driven by combined direct and alter- 
nating currents I cx t{t) = Id + Ia^iiiLut, fractional giant 
Shapiro steps at voltages 



(V) 



2eJ\q 



(1) 



where (■ ■ •) is the time average and n is an integer, have 
been observed both in experimentsa and in computer sim- 
ulationsia Qualitative arguments have been proposed 
to esglain these fractional steps iu,i.erms of vortex mo- 
tioroQ and topological invariancefm In addition to the 
fractional steps, a scries of small subharmonic steps has 
been found .E3 

Two slightly different models have been used to catch 
the essential properties of a JJA: the resistively shunted 
junction (RSJ) model aiad the time-dependent Ginzburg- 
Landau (TDGL) model.0 The RSJ model is based on the 
assumption that all the current goes through the array 
and that the current is conserved locally at each instant. 
The TDGL model in the absence of an external current 
describes either a situation where all the current goes 
through the array, but where the local current conserva- 
tion is relaxed, or a situation where not all the current 
goes through the array (leakage to the ground) and the 
current is conserved at each instant. The former view 
means that the TDGL model can be regarded as a sim- 
plified version of the RSJ model and at the same time as 
a less restrictive model of a JJA. The latter view has led 



to the suggestion that a JJA with localjiamping is a pos- 
sible realization of the TDGL modelJijtj In the presence 
of an external current the physics of the TDGL model 
depends on the choice of the boundary condition. Wc 
use here a boundary condition corresponding to the case 
when the normal current flow is through the array just 
as in the RSJ case. However, current is only conserved 
globally and not locally. In this way, we can compare the 
effects caused by the difference in local current conserved 
dynamics as in the RSJ case with the TDGL dynamics 
which only has global current conservation.E.3 

Both dynamic models are equivalent as far as static 
equilibrium properties are concerned, sinae-jthey have the 
same equilibrium Boltzmann distribution.E£l On the other 
hand, for dyrpxuic quantities such as the-dynamic dielec- 
tric functionJl3'E£l flux-noise s^p^jdjrumJi^O and current- 
voltage {I-V) characteristics,t§t3 the equivalence is not 
guaranteed. It has recently been suggested that the 
TDGL model could describe the fhpi-noise experiment 
for a JJA better than the RSJ modelJl3 However, a some- 
what different conclusion was reached in Ref. [15| where 
properties like the linear response and nonlinear I-V 
characteristics were found to be the same for the two 
models. 

In Ref. [l5] a novel boundary condition [the fluctuat- 
ing twist boundary condition (FTBC). based on global 
current conservation was introduced.EJa We show in this 
work that the very same Shapiro steps are found in the 
TDGL and RSJ models .when we employ the FTBC as 
the boundary condition. til This suggests that the exis- 
tence of the steps does not depend on the details of the 
dynamic models: This robustness can be explained by 
the topological nature of steps where the ground state 
degeneracy (both models are equivalent in_this respect) 
has been shown to play an important roleJj The widths 
of Shapiro steps as a function of I a and u have also been 
a subject of much interest. For example, the maximum 



width of the integer and the fractional steps have been 
shown to have-a-different frequency dependence for larger 
frequencies.B : E3fcy We find that the maximum width of 
the half-integer step for the TDGL dynamics has a dif- 
ferent frequency dependence than for the RSJ case. This 
offers an experimental possibility (similar to the experi- 
ment on an / = array in Ref. B2h to investigate whether 
a JJA could sometimes be better described by the TDGL 
dynamics. 

In the presence of applied direct currents, the critical 
currents I c (f), beyond which the voltage takes nonzero 
values, rbaye been measured in experiments! 2 . 3 ! arid simu- 
lations.QCfJ-Although theoretical predictionsBH and ex- 
perimentsE 3 ! for / = 1/2 give the result I c (f — 1/2) = 
\/2 — 1 ps 0.414 in units of I c (f = 0), computer simu- 
lations with the conventional mcthcd_nf_nniform current 
injection gives the value 0.35 ± 0.01.1 LlEfl In Ref. | it has 
been argued that this discrepancy is due to the bound- 
ary condition and that the conventional method destroys 
the translational symmetry of the ground state. On the 
other hand, we find in this work that the FTBC gives 
the value 0.4142(1) for both the RSJ and TDGL mod- 
els, which suggests that the FTBC is a more adequate 
boundary condition since it conserves translational sym- 
metry. 

We start by introducing the equations of motion for the 
RSJ and TDGL models with the FTBC (see Ref. |l| for 
details). In the FTBC the twist variable A = (A x ,A y ) 
is introducecLand the gauge-invariant phase difference is 
changed intccJ 

fa =9i- 9j - Aij - Tij - A, (2) 

where 9i is the phase of the superconducting order pa- 
rameter at site i, ry = r; — tj is a unit vector from site 

i to j, and Ay = (2-7r/$o) /• A • dl with the magnetic 
vector potential A and the flux quantum <E>o for Cooper 
pairs. 

In the RSJ model, the equations of motion for phase 
variables are determined by the local current conservation 
at each site (see, for example, Refs. [| and |ic| ) : 

0i = - ^2 G v^2 sin ^ ( 3 ) 

j k 

where the primed summation is over four nearest neigh- 
bors of j, Gij is the square lattice Green function, and 
the unit of time is h/2eRI c with shunt resistance R and 
critical current I c of the single junction. In this work 
we only consider the array at zero temperature and ac- 
cordingly the thermal noise terms are disregarded (see 
Ref. |5| for finite temperatures). Foii_tlip TDGL model 
the equations of motion are given byEjIlj 

0i = sin^y, (4) 

3 

where t is in units of 2e/I c . 



In the FTBC case, the periodicities of the phase vari- 
ables are preserved in both directions, i.e., 9i = #i+L X = 
®i+Ly, and thus the voltage drop in each direction across 
the whole array is given by 

Ri. HL ■ 

from the Josephson relation. From the condition of global 
current conservation in each direction, we obtain the 
equations of motion for the twist variables of the array 
driven by axL.external current Text (in units of I c ) in the 
x direct ion :Ej 

dA x 1 ^-^ 

= — df + L2 ^ sln Wi> ( 6 ) 

dA v 1 \ , , , 

where S(tj>« (E(ij> ) i s the summation over all nearest- 
neighboring pairs in the x(y) direction. 

It should be observed that we are here considering the 
TDGL model [as specified by the dynamical equation (0)] 
with globally conserved current [as specified by Eqs. (|6|) 
and ([fl)]. One can also consider the TDGL model [as 
specified by Eq. (|J) ] without global current conservation 
within the plane [i.e., without Eqs. @ and (0)] aneLin- 
terpret it as a model with leakage to the ground.l!3ll 3 J In 
this case an applied external current only leads to dissi- 
pation to the ground at the boundaries where the current 
is injected and extracted. Consequently there exist no gi- 
ant Shapiro steps for the TDGL model without globally 
conserved current within the plane Ej 

We first consider an array with an external current 
-^cxt = Id + la siriLut. We use the Euler algorithm with 
discrete time step At = 0.05 to integrate the equations 
of motion [Eqs. (0), (^|), and (^) for the RSJ model, and 
Eqs. (|), (§), and (0) for the TDGL model], and the time- 
averaged voltages (V) in units of LUu>/2e are calculated 
from Eq. (||). We adopt the simulated annealing Monte 
Carlo method to find the ground states of the array and 
then use them as initial conditions of phase variables to- 
gether with A x (t = 0) = Aj,(0) = 0. Figure shows 
the fractional giant Shapiro steps in the I-V character- 
istics for (a) the RSJ model and (b) the TDGL model. 
Although we have found quantitative differences in the 
Shapiro steps for the TDGL model (e.g., small step sizes 
at (V) = 1/3 and 1/2 for / = 1/3 and 1/4, respectively), 
it is clear that with the FTBC not only the RSJ but also 
the TDGL model generates the integer and fractional 
steps. We have also observed weak subharmonic steps 
for both models as in Ref. [h]. 

Figure [2] displays the maximum widthsEJt3 of the steps 
(V) = 1/2 and 1 versus u> for both models. Although 
both show the same qualitative behavior in the low- 
frequency regime, it is apparent from the figure that 
the high-frequency behaviors of the half-integer steps are 
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different for the two models. Since the frequency de- 
pendence of the maximum width can be measured for a 
J JA,E3 this offers the possibility of experimentally distin- 
guishing between the two types of dynamics. 

We have also performed computer simulations apply- 
ing a constant direct current J ex t (t) = Id for three differ- 
ent cases and obtained the critical currents I c (f)' One 
case is the RSJ model with the conventional method 
of uniform current injection which employs the periodic 
boundary condition (the free boundary condition) in the 
direction perpendicular (parallel) to the applied currents. 
The other two are the RSJ and TDGL models with the 
FTBC. We present in Table | a comparison of these three 
cases, which reveals that the FTBC gives correct val- 
ues for both the RSJ and TDGL models. For the con- 
ventional current injection method we obtained different 
values, e.g., I c (f = 1/2) = 0.35(1), as was also found 
in Ref. [|. We checked the system size dependence for 
L = 4, 8, • ■ ■ , 128 and found no change. Nevertheless, 
these smaller values are caused by the boundary condi- 
tion which destroys the translational symmetry of the 
ground state.u In Ref. |^ this problem was circumvented 
by a nonuniform injection method which matched the 
translational symmetry of the ground state and the cor- 
rect value 0.414 was found. In the FTBC case the trans- 
lational symmetry of the ground state is automatically 
preserved and consequently this boundary condition di- 
rectly yields the correct result. We also-calculated the 
critical current with the busbar geometrytS and obtained 
an even smaller value of I c (f = 1 /2) than for the conven- 
tional uniform injection method, as was already noticed 
in Refs. |5| and |24|. From these comparisons we conclude 
that the FTBC has an advantage over other commonly 
used boundary conditions. 

In Fig. H we show the average energy defined by 
E = ~(J2{ij) cos (j)ij (t)) / L 2 as a function of Id for sev- 
eral cases. Our results for the RSJ model with the FTBC 
are in perfect agreement with the results from the ana- 
lytic equations given in Ref. ^0[ which suggests that our 
results contain no boundary or finite-size effects. The 
TDGL model is found to give the same E for currents less 
than the critical value. Beyond the critical current, the 
TDGL model gives a lower energy, implying that the ar- 
ray is closer to the ground state than the RSJ model. We 
believe that this explains the robustness of the 1/2 step 
of the TDGL model at high frequency (see Fig. g), since 
it is expected that the ground state and its vortex super- 
lattice structure ijdays an important role in creating the 
half-integer steps.H We find in all cases that E(Id) has a 
cusp structure at the critical current. One may also note 
in Fig. U that the conventional uniform current injection 
method leads to a result which differs from the exact an- 
alytical result. Figure |] gives the critical currents I c (f) 
at / = p/q with q = 1, 2, • • • , 8 (for comparisons with pre- 
vious works, see Ref. |^). For all values of /, we obtain 
identical values of I c (f) for the TDGL and RSJ models. 

In conclusion, we have performed simulations for the 



RSJ and TDGL models subject to the FTBC. Fractional 
giant Shapiro steps are obtained for both models, which 
suggests that the existence of the steps does not depend 
crucially on the condition of instantaneous local current 
conservation. However, the maximum width of the half- 
integer step at / = 1/2 has a qualitatively different high- 
frequency behavior for the two models. The critical cur- 
rents of the array with direct applied currents were also 
calculated for both the models subject to the FTBC and 
compared with the results obtained for the RSJ model 
with the conventional method of uniform current injec- 
tion. It was concluded that the FTBC for both models 
gives values in agreement with experiments and analytic 
calculations, while the conventional method fails in this 
respect. 

The present calculation supports the conclusion 
reached in Ref. [l^ that the TDGL and RSJ models with 
the FTBC are qualitatively equivalent for low frequen- 
cies (compare Fig. ^|) and small currents (compare Fig. |^) 
whereas for larger frequencies and larger currents there 
exist qualitative differences. The fact that both the mod- 
els have qualitatively similar sequences of giant Shapiro 
steps suggests that the existence of these steps is strongly 
linked to an. equilibrium property like the ground state 
degeneracy!] 
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TABLE I. Comparison of critical currents at / = 0, 1/2, 
and 1/3 for the RSJ model with the conventional method 
(7™ nv ), the RSJ model with the fluctuating twist boundary 
condition (1? SJ ), the TDGL with the FTBC (/J DGL ), and 
the analytic results in Ref. ^ (7c nal ). All values are in units 
of I c (f = 0) and the numbers in parentheses are numerical 
errors in the last digits. 



/ 


rconv 


rRSJ 


rTDGL ranal 





1.00(1) 


1.0000(1) 


1.0000(1) 1 


1/2 


0.35(1) 


0.4142(1) 


0.4142(1) 0.41421 


1/3 


0.14(1) 


0.2679(1) 


0.2679(1) 0.26789 
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Id 



FIG. 1. Time- averaged voltages (V) in units of Lhui/2e 
versus direct current Id for (a) the RSJ and (b) the TDGL 
models in an L x L Josephson junction array in case of / = 0, 
1/2, 1/3, 1/4, 1/5, and 2/5 (from the left to the right). The 
sizes of arrays are L — 4 (for / = 0), 8 (/ = 1/2 and 1/4), 
9 (/ = 1/3), and 10 (/ = 1/5 and 2/5), and we have used 
the fluctuating twist boundary condition together with the 
condition of global current conservation and applied external 
currents J cxt = Id + la sin cut with I a = I c and u>/2n = 0.1 in 
units of 2eRI c /h for the RSJ model and I c /2e for the TDGL 
model, respectively. Fractional giant Shapiro steps are clearly 
shown for the TDGL as well as for the RSJ model. All curves 
except / = are horizontally displaced for clarity. 



(V)= 1/2 (TDGL) 




FIG. 2. Frequency dependence of the maximum widths of 
the Shapiro steps at (V) = 1/2 and 1 for the RSJ and TDGL 
models with the FTBC in case of / = 1/2. Our results for the 
RSJ model are in good agreement with the analytic results. 
(Ref. ^). The high-frequency behavior of the 1/2 step for 
the TDGL model is shown to differ from the RSJ model. The 
lines are guides to the eye. 
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FIG. 3. Energy per site, E = — (Xwy) cos <f>ij)/L 2 , for an 
L x L array in case of / = 1/2 with an applied direct current 
Id- The results from the analytic equations (solid curve) in 
Ref. |2(] are in perfect agreement with our results for the RSJ 
model using the FTBC. The TDGL model gives the same E 
below the critical current. However, the RSJ model with the 
conventional method has a cusp structure at a different value 
of J*. 
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FIG. 4. Critical currents I c (f) as a function of / for the 
RSJ (denoted by +) and the TDGL (o) at / = p/q with 
q — 1, 2, • • • , 8. The fluctuating twist boundary condition is 
used together with the condition of global current conserva- 
tion. For all values of / tested in this work, the RSJ and 
TDGL models give the same value of I c (f) within numerical 
accuracy. The line is a guide to the eye. 



G 



